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Abstract. Constraints on the density and thermal 3D structure of the dense molecular cloud core p Oph D 
are derived from a detailed 3D radiative transfer modeling. Two ISOCAM images at 7 and 15 /im are fitted 
simultaneously by representing the dust distribution in the core with a series of 3D Gaussian density profiles. 
Size, total density, and position of the Gaussians are optimized by simulated annealing to obtain a 2D column 
density map. The projected core density has a complex elongated pattern with two peaks. We propose a new 
• method to calculate an approximate temperature in an externally illuminated complex 3D structure from a mean 

optical depth. This "T— "-method is applied to a 1.3 mm map obtained with the IRAM 30m telescope to find the 
approximate 3D density and temperature distribution of the core p Oph D. The spatial 3D distribution deviates 
strongly from spherical symmetry. The elongated structure is in general agreement with recent gravo-turbulent 
collapse calculations for molecular clouds. We discuss possible ambiguities of the background determination pro- 
cedure, errors of the maps, the accuracy of the T^method, and the influence of the assumed dust particle sizes 
and properties. 
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1. Introduction 1.1. Verification of ID models of early star formation 

by observations of molecular cloud cores 

Stars and planetary systems are believed to form within . . ... 

, ,, , , , , . ,, r j In the case ol low-mass star-forming regions, it has tra- 

mterstellar molecular clouds, irom the collapse ol conden- „„'... 

„ „, „ , . „ ditionally been thought that stars form from the mside- 

sations ol matter often called prc-stellar cores. Despite „ n . ° . . . . 

, , . , , ,, , . , , . , out collapse oi singular isothermal spheres that are mi- 

rcccnt observational and theoretical progress, the physical .... . ? .... . i M1 1 — ri| ^o4 , 

i j. , ,i r • tii. j ,i • tially m quasi-static equilibrium bhu ct al. 1987 , sup- 

processes leading to the lormation ol these cores, and their , r — \ — \ — — \ — 

„ 11/ H — — , 1 1 loorinl\ ported against gravity by magnetic and thermal pressure 

collapse remain poorly known (e.g. Andre ct al. 2000). . , , . ° . _ . 

, , „ , „ . . , , , and evolve only d ue to slow ambipolar diilusion processes 

I he study of pre-stellar cores is essential to the under- r— ■ : — i t' " i , 

£ , r .. . Mouschovias 1991). In this case, prc-stcllar cores at the 

standing ol star lormation, since these objects represent * — — " 7™ — * , „ „ . , . 

r i ii i • i • a onset ol collapse arc expected to follow a l/r z density 

the initial conditions ot gravitational collapse which mfiu- „ ... i " r i r n „ i , ' . , 

, r , , r prohle at all radii r bhu 1977). I hen they collapse dy- 

ence the subsequent stages ol the star-lormmg process. 1 . . — 7, , . 

namically under their own weight and form a protostar, 

which subsequently accretes matter from its surrounding 

envelope. 

In the last decade, much effort has been put in the un- 

Send offprint requests to: stein@mpia.de derstanding of these cold (T ~ 5-15 K) pre-stellar objects, 
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which emit most of their radiation in the (sub)millimeter 
domain with observations that benefitted from the de- 
velopment of millimeter bolometers, sensitive line de- 
tector systems, and infrared cameras. Those observa- 
tions made it possible to derive line-of-sight integrated 
properties of the pre-stellar cores. Several studies using 
(sub-)millim e ter continuum llWard-Thompson et alll994t 
Andre et all Il996t iMotte et all Il998l lLaunhardt et al 



2004) or mid-infrared (M IR) absorption l|Bacmann et al 



200rl lAbergel et alllflflrfl have in particular addressed the 



density structure. Determining the density structure of 
pre-stellar objects is fundamental to understand core and 
protostar evolution, since it has been previously shown 
that the density structure of pre-stellar objects influences 
the proto-stellar accretion rat e, i.e. after the formatio n of a 
cent ral object. According to Henriksen et al.1 l)l997(l (see 



rdmg 

also iFoster fc Chevalier! 1 1993() . non-singular density pro- 
files in pre-stellar cores lead to accretion rates declining 
with time in the proto-stellar stage, whereas the accretion 
rate is e xpected to be constant if th e core has singular r~ 2 
profiles llShulll977t Ishu et al.lll987h . 

It was found from the mm and MIR maps that column 
density profiles of pre-stellar objects are flattened out in 
their centers. Interpreted by ID singular isothermal sphere 
models, the derived density distributio ns do not follow the 
r _2 -profile from the IShu et all l)l987(l model, though the 
profiles do have an r~ 2 part. 

The fact that some co res even present sharp edges at 
large radii (~ 0.1 pc, see iBacmann et al.ll2000|) tends to 
show that the reservoir of mass available for subsequent 
accretion is limited. Therefore to some extent, the density 
structure of pre-stellar cores must gover n stellar masses 
and influence the initial mass function l|Umemoto et all 
120031) . 

Using the ID approximation, several types of hydro- 
static and quasi-static models have been proposed to ac- 
count for core formation and evolution. Comparison of 
the density profiles derived from the observations with 
the models has been used in order to discriminate be - 
tween various models (e.g. IWa rd-Thompson et al.lll99 
Ward-ThompsorJll998tlAndre et al.lll996HBacmann et al 



2000). In most of these studies, circular (or elliptical) av 
erages of the dense core emission (either in the mm or in 
the MIR) have been performed according to the apparent 
(projected) geometry of the object. The azimuthal aver- 
age was also used to increase the signal-to-noise ratio. The 
profiles thus derived were compared to various power law 
density models. It was found that the ambipolar diffusion 
models of magnetically-supported cores provide the best 
fits, requiring a background m agnetic field ranging from 
30 to 100 fiG llWolf et al.ll200l . 



1.2. Limitations of ID models of early star formation 

Both from observational and from the theoretical side, 
there are many indications that star formation is an in- 
trinsically three-dimensional process. The profile analysis 



was often performed in ID or 2D because of the difficulty 
of obtaining information along the 3rd dimension, as well 
as the substantially enhanced numerical effort in 3D star- 
formation simulations and 3D radiative transfer modeling. 

However, the increased sensitivity of bolometer ar- 
rays as well as the observations of cores in absorption 
in the mid-infrared (a technique whic h is more sen- 
sitive to small column den sities - e.g. IBacmann et alJ 
|200(A ISteinacker et~ai1 12004|) reveals filamentary struc- 
tures for the cores at large scales, which significantl y 
differ from near-spherical approximations I Aoail l2004h . 
Stati s tical analyses of core projected shapes ijjones et all 
12001b |j ones also tend to show that cores 

are triaxial. If ID interpretations represent in simple cases 
(i.e. for isolated cores) a good approximation of the cores, 
they can a lso lead to erroneous densi t y stru ctures, as illus- 
trated bv_ B^ligg^g^^Pa^^g^^^ l|2003j) . More specif- 
ically, ISteinacker et al.1 l|2004h used non-symmetric, fila- 
mentary density distributions from a smooth particle hy- 
drodynamics simulation of a cloud core evolution to pro- 
duce MIR and mm maps. By interpreting the absorption 
maxima in MIR maps and emission maxima mm maps 
with a ID model, they found density profiles in ostensi- 
ble agreement with standard ID core formation models, 
but in contradiction to the true underlying non-symmetric 
distribution. 

Moreover, the increased sensitivity and resolution of 
upcoming telescopes will make it possible to see small- 
scale structures in these objects and make more sophisti- 
cated models more and more necessary. The 3D density 
structure and eventually its evolution in time is also of cru - 
cial interest for chemical models (e.g. lAikawa et ai1E0 03). 



1.3. Gravo-turbulent star formation 

The hypothesis that star formation mostly occurs in 
magnetically subcritical cores and is regulated by slow 
ambipolar diffusion processes llShu fc Adams! Il987l) has 
bee n challenged f rom several sides (see, e.g. , the reviews 
by iLarsenl 12003 iMac Low & Klessenl \2004 . It is only 
applicable to isolated, single stars, while it is known 
that the majority of stars form in small aggregates 
or la rge clusters llAdams fc Mversl l200lt lLada k, Ladal 
2003). Furthermore, t here is both observational evi- 



l200ll iHenning et al.ll200ll IWolf et alJ 


2003) and theoret- 


ical reasoning (e.g. iNakand 


1998 


that cloud cores do 



not have magnetic fields strong enough to support the 
core against gravitational collapse. Also, the long life- 
times implied by the quasi-static phase of evolution in 
the model are currently discussed r egarding, e.g., ob- 
servational statistics of cloud cores l)Tavlor et all Il996t 



iLee & Mversl 


1999: IVisser et all 12002 


) and chemical age 


considerations (Ivan Dishoeck & Blake 


Il998l lLanger et alJ 



200(J. iTassis fc Mouschoviasl l(2004h . however, have ar- 
gued that the observational data are consistent with 
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the predictions of the theory of self-initiated, ambipolar- 
diffusion-controlled star formation. 

An alternative to mediation by magnetic fields as con- 
trolling age nt for stellar birth is sup ersonic interstellar 
turbulence ( Mac Low &: Klessenl 1200^1 . The presence of 



highly supersonic (and usually super- Alfvenic) turbulent 
motions in Galactic gas clouds is inferred from molec- 
ular line measurements. The observed line widths typi- 
cally exceed__the values for thermal line broadening by 
far l|Blitzlll993ft . Only on very small scales and for the 
so-called quiescent cores does the total linewidth become 
comparable to the thermal width alone. It is important 
to note that turbulence usually carries sufficient energy 
to counterbalance gravity on global scales. The presence 
of supersonic turbulence therefore prevents contraction of 
the cloud as a whole. On small scales, however, it may 



ElmcE'rccn 


1993; Padoan 1995; Ballcstcros-Paredcs ct al. 


1999; Klcssen et al. 2000; Padoan & Nordlund 1999, 2002; 


Heitsch et alJ 2001 


). Supersonic turbulence establishes a 



complex network of interacting shocks, where converging 
flows generate regions of enhanced density. The compres- 
sion can be sufficient for local gravitational instability to 
set in. The same random flow that creates density en- 
hancements in the first place, however, can also disperse 
them again. For star formation to occur, the localized com- 
pression must be fast enough for the contraction region to 
'decouple' from the flow. 

The theory of gravo-turbulent star formation pre- 
dicts that molecular cloud cores form at the stag- 
nation points of the complex turbulent flow pattern. 
Their density profiles are constant in the inner region, 
exhibit an approximate power-law fall-off further out, 
and usu ally appear truncated at som e maximum ra- 
dius fsee lBallesteros-Paredes et al.ll2003l) . Such configura- 
tions c l osely resemble equil ibrium Bonnor-Ebert spheres 
l)Eherd 1-195.4 iBonuo^ ITflSfih . which arc popular models 
for interpreting observed cloud cores. The best studied 
example probably is t he Bok globule B68 ijAlves et alJ 
l200lULada et al.ll2003h . However, supersonic turbulence 
does not create hydrostatic equilibrium configurations. 
Instead, the density structure is transient and dynami- 
cally evolving, as the diff erent contributions to virial equi - 
librium do not balance ijVazouez-Semadeni et alJl2003|) . 
Numerical models predict that gravo-turbulent fragmen- 
tation of molecular cloud material will lead to cores that 
are generally triaxial and in many cases highly distorted. 
Depending on the projection angle, they often appear ex- 
tremely elongated, being part of a filamentary struct ure 
that may connect several objects l|Klessen et al-lEoO^I . 

It is the goal of this paper to extend our understanding 
of pre-stellar cores by entering the complex and difficult 
world of 3D modeling. The target of our investigation is 
the pre-stellar core p Oph D located in the p Ophiuchi 
star-forming cluster. This core was selected because it has 
been observed both in the MIR and mm range and the 
maps indicate a cloud structure deviating from spherical 
symmetry. We start by outlining the general strategy of 



the analysis in Sect. 2. In Sect. 3, we present results of 
a simultaneous fit of ISOCAM images at 7 and 15 pm 
by a series of 2D Gaussian profiles to obtain the number 
column density at each image pixel. The necessary back- 
ground determination scheme is described and discussed. 
We compare the derived dust density to an IRAM 30m 
1.3 mm map in Sect. 4. Presenting and applying a new 
method to obtain 3D density and temperature informa- 
tion of cores, we derive constraints on the structure of the 
core p Oph D. A discussion of ambiguities and errors of 
the obtained distributions is given in Sect. 5. In Sect. 6, 
we compare the derived spatial distribution with distri- 
butions typically obtained in gravo-turbulent simulations 
and conclude our findings in Sect. 7. 



2. General strategy of the 3D analysis 

As the 3D analysis of molecular cloud cores based on 
continuum images is new and includes several steps and 
approximations, we present in this section a simplified 
overview of the applied method. Details, tests, and fur- 
ther discussion of single points have been moved to the 
later sections for clarity. 

First, we determine the column density of the core re- 
gion for each pixel of the 7 and 15 /xm images. In order 
to estimate the background radiation that is absorbed by 
the core, we consider the un-absorbed background radia- 
tion near the core and interpolate it behind the core. The 
optical depth is calculated for each pixel and, under the 
assumption of dust particle size and opacities, converted 
to a dust column density. To model the 3D distribution of 
the core, we use a series of 3D Gaussian distributions with 
axes aligned to the Cartesian axes. The position, the half 
width at half maximum (HWHM) along the three axis, 
and the density weight contain 7 free parameters for each 
Gaussian. To model the column density, the position in 
the plane of sky, two HWHM, and the weight are varied 
to fit the images. The corresponding x 2 -function is mini- 
mized using simulated annealing. 

The dust density and temperature distribution along 
the line of sight is determined by considering mm emission 
of the dust particles absorbing in the MIR. The emission 
seen in the 1.3 mm image is obtained from a line-of-sight 
integral over the emission at each location in the core de- 
pending on both the local dust density and temperature. 
By varying the position and HWHM of the Gaussians 
along the line of sight, we can fit the 1.3 mm image to 
obtain the 3D density and temperature structure apply- 
ing again the simulated annealing optimization algorithm. 

The calculation of the dust temperature distribution 
for a given density distribution requires the application of 
a 3D continuum radiative transfer (CRT) code. As these 
computations are far too time-consuming to be repeated 
many times, we propose and apply an approximative re- 
lation to calculate the temperature T at a location in the 
core from a mean optical depth f at 7 pm obtained by a 
weighted averaging over all directions. This T^relation is 
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obtained from a few 3D CRT runs covering average and 
limiting density distributions. 



3. Determination of the dust number column 
density map from ISOCAM maps at 7 and 15 

3.1. Background determination 

To obtain a detailed dust column d ensity map of the core , 



we use the ISOCAM maps from Abcrgcl ct al. 
which also have been discussed in lAbergel et al 



1996), 



1998) 



The maps were observed at a central wavelength of ~ 7 /im 
(LW2 filter [5-8.5/tm]) and 15/im (LW3 filter [12-18 /an]) 
with a resolution of 6". The MIR emission seen in the p 
Oph cloud is thought to arise from very small grains or 
Polycyclic Aromatic Hydrocarbons (hereafter PAHs) ex- 
cited on the outside of the cloud by far-UV radiation. The 
PAHs then re-emit the absorbed energy in the form of 
spectral bands in the MIR, of which several are encom- 
passed in the ISOCAM filters. Since the extinction in the 
far-UV is very high and does not penetrate further than 
Ay ~ 1 mag, the core itself should be free of PAH emis- 
sion, the bulk of the observed MIR emission arising from 
the outer parts of the cloud. For the same reason, emission 
of the cold dust particles within the core in the 7- 15/im 
domain is highly unlikely (the core temperature is lower 
than 20 K). The radiation received can be modeled as fol- 
lows (see also Fig.0: 1. A background layer at the origin 
of excited PAH emission from the outer cloud parts (la- 
beled bg), 2. The dense core absorbing this background 
radiation in the MIR (labeled c), 3. A foreground layer 
similar to the background layer containing absorbing cold 
dust and a layer of emitting small particles at the border 
of the cloud (fg), and 4. A layer of zodiacal dust with 
emission and absorption (zd). With this approximative 
decomposition, the total intensity along the line of sight 
can be written as 



obs 



,-Tc „ — 



I 



fg 



(1) 



where r is the optical depth. 

The p Oph cloud i s illuminated from behind by 
the B2V-star HD147889 l|Liseau et al.lll999|) exciting the 
PAHs of the background layer. As there is no star in the 
vicinity of the cloud illuminating the foreground layers, 
the PAH emission arising from the foreground layer is sup- 
posed to be negligible. This is in agreement with the re- 
sults of Bacmann et al. (2000) who found that the entire 
foreground emission came from zodi acal dust. The zodia- 
cal dust emission was determined bv lAbereel et alj l)l996|) 
to be 5.7 and 39 MJy/sr for the LW2 (7/zm) and LW3 
(15 /im) maps, respectively, and has been subtracted from 
the map. Currently, a possible absorption by cloud dust in 
front of the core can not be distinguished unambiguously 
from the core absorption. We note that a small part of the 
dust distribution derived from the maps might belong to 
the foreground cloud rather than to the core itself. In that 



c fg zd 





Fig. 1. Sketch of the possible layers of matter altering the 
background radiation on its way through the dense core. 



sense, the derived dust densities will be an upper limit to 
the real densities. 

To determine r c and the dust column density map, the 
background behind the core has to be interpolated from 
the observed background in the vicinity of the core. We 
assume that the background behind the core follows the 
gradients that are visible in the maps. After excluding 
bad pixels containing contamination from previous obser- 
vations and periodic remnants, we determine a first back- 
ground approximation by manually masking the core area 
and interpolating the outer background into the masked 
area. The background at a given pixel was interpolated 
as the average of all non-masked pixels within a given ra- 
dius. This radius was chosen to be 1/3 of the map extent. 
This reproduces well the background gradient throughout 
the map without emphasizing local noise features. After 
having defined this first-order background, we redefine the 
core area to be the range of pixels where the absorption 
of this background is larger than 1/e. Again, the nearby 
background is interpolated into this new core area. We 
have verified the stability of the results of our fitting pro- 
cedure by varying the averaging radius and the size and 
shape of the initial manually-determined mask of the core 
area. 

3.2. Determination of the optical depth and column 
density map 

The interpretation of the obtained optical depth along 
lines of sight through the core is ambiguous. The absorp- 
tion coefficient of dust particles K is a function of their 
chemical composition C, of their size a, the temperature 
T, the wavelength A, and (via C, a and T) the location 
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in space. As th e temperature range of roughly 5-20 K is 
small (see, e.g.. IStamatellos et al.ll200 l), we neglect tem- 
perature dependencies for the absorption coefficient. As 
far as dust size is concerned, the dust may grow in the 
innermost parts of the core changing both its size and 
chemical composition and hence its optical properties. A 
mixing of the dust particles due to gas motion countering 
spatial dependencies is possible but hard to model due 
to the complex turbulent motions in these cores. A mini- 
mum modeling of the dust evolution and the gas m otion 
would be required (see lOssenkopf k. Hennind 1994), but 
it would introduce a wealth of more parameters into the 
fitting procedure. Since in this paper, we would like to em- 
phasize the detailed structure that was often neglected in 
former models, we assume that absorption is dominated 
by dust of a typical size a and refer to a forth-coming pa- 
per dealing with a possible spatial variation of the dust 
size and absorption coefficient. 

The dust number column density is calculated by ray- 
tracing the background intensity through a given spatial 
3D dust number density distribution representing the core. 
There are several useful basis functions that can be used to 
describe the spatial 3D dust density distribution n(x, y, z) 
of the core in a series expansion in Cartesian coordinates. 
Here, we choose a series of g 3D Gaussians with their main 
axis aligned with the coordinate axis, having the form 



n(x, y,z) = ^2 G i ex P [- d i( x > z )\ 

i=l 

with the argument 



(2) 



di{x,y, z) 



(3) 



The parameterization contains 7 * g free parameters 
(weights Gi, translations Xi, yi, Zi, and size parameters 
&i, bi, c{). This representation has the advantage that the 
Gaussian profiles are a good approximation for the flat- 
tened innermost part of t he cores that is seen for many 
cores (for profiles see, e.g.. iBacmann et al.l EoOO'l . and in- 
clude the spherically symmetric case. For the outer parts, 
the gradients of the density profile have a large uncer- 
tainty as projection effects severly change the profiles that 
are obtained from co lumn density modeling (see, e.g., 
ISteinacker et al.l [2004) . Therefore, it is currently uncer- 
tain if the density follows a powerlaw or an exponential. To 
keep the inversion method proposed in this paper numer- 
ically feasible, we have chosen the number of Gaussians 
g to be 30, and the implications of this choice are dis- 
cussed in Sect. 4.2 and Sect. 5.1. We also have tested 
a series of Gaussians with axes inclined to the line of 
sight and got the same overall column density distribution 
with a lower number of Gaussians. The numerical effort 
is higher, however, since the formula for Gaussians with 
axes inclined to the line of sight is substantially more com- 
plicated. Moreover, for a given number density n{x,y,z), 



conveniently, the number column density 

N {y,z) = dx n(x,y,z) 



Vtt^ a-iGi exp [-d(xi, y, z)} 



(4) 



i=l 



can be calculated analytically if the line of sight is cho- 
sen to be, e.g., along the x-axis. An interclump medium 
represented by a constant density component has not been 
added to the density as the Gaussians include this case for 
large size parameters. Assuming dust particles with an in- 
ternal density of pd, the corresponding total dust mass in 
the core is 



M = -TT 5/2 a 3 p d Gja^iCj. 



(5) 



The optical depth between two points x\ and x-i on the 
line of sight is 



T x (xx,x 2 ,y,z, A) = <r(A) dx n{x,y,z) 



(G) 



with the absorption cross section <r(A). For a ray crossing 
the entire core, we find 



t x (-oo, oo, y, z, A) = a(X)N(y, z). 



(7) 



Summing over all pixels k and maps 1, we define the x 2 - 
function of the fit to be 



m Ni 
1 = 1 k=l 



(8) 



where Di contains the flux values of each pixel of the two 
images, Fi is the background flux after undergoing absorp- 
tion due to the model core, m — 2 is the number of maps, 
and Ni is the number of pixels of map /. 

The optimization of the 150 parameter dust number 
density by fitting about 1430 pixels of the MIR maps was 



algorithm fairkpatrick et alJ 


Il983t IVanderbilt & Louiel 


1984; 


Thamm et alJ 


1994; 


Steinacker & Thamm 


1996; 


Steinacker & Hennine 


II2003I) which is able to both find 



minima in the manifold parameter space and to find the 
deepest minima. In Fig. [5J the relative location and size 
of the individual Gaussians is indicated by ellipses having 
HWHM of size di and bi . The line thickness of the ellipses 
is logarithmically proportional to the density weight of the 
Gaussians. The general form of the core is an elongated 
structure with two peaks and a ratio of projected axis of 
about 4:1. The larger southern peak is well described by 
three almost spherical Gaussians. The northern peak has 
a complex shape and a correct representation of its den- 
sity structure may require more Gaussians than used in 
this analysis. 

The number column density map is shown in Fig.|21as a 
function of the pixel index of the map (one pixel has a size 
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10 20 30 40 50 

w pixel 



Fig. 2. Relative locations and sizes of the individual 
Gaussians for the optimized number column density dis- 
tribution. The ellipses are cuts of the Gaussians for x = Xi 
and rii — Gi/e 2 . The thick-lined ellipses indicate the most 
dense Gaussians. 



of 960 AU). The double-peaked structure is pronounced, 
with the southern peak being the denser part. 

In Fig. the optical depth map for 7 /im is shown 
as a function of the pixel number. The contours are la- 
beled with the optical depth value t x along rays com- 
pletely crossing the core, reaching a maximum of about 
0.5 in the southern maximum and less than 0.4 in the 
northern maximum. The optical depth does not exceed 
0.14 in the region between the two maxima. 

The resulting fit maps using the optimized number 
density distribution are given in Fig. [5] along with two of 
the original ISOCAM images (upper and lower left panel 
shows 7 and 15 /im, respectively). The white line in the 
upper left panel represents a spatial scale of about 0.13 
arcmin or 0.1 pc (20000 AU). The fit is not able to repro- 
duce the small-scale noise that is part of the real maps but 
resembles well the overall appearance of the cloud core. 

With 30 Gaussian basis functions involved, the fit is 
ambiguous. The most obvious ambiguity when running 
the optimization procedure arises from identical density 
structures where Gaussians just have been exchanged. 
Moreover, there are several possibilities to arrange the 
Gaussians to represent the complex extended double-peak 
structure of such a core. These correspond to different lo- 
cal minima in the parameter space and it requires a special 
minimization scheme like the simulated annealing algo- 
rithm as used in this modeling to find the global minimum 




10 20 30 40 

§ pixel 

Fig. 3. Map of the optical depth t x at 7 /j,m of the cloud 
core derived from simultaneously fitting ISOCAM maps 
at wavelengths of 7 and 15 /im. Each pixel has a size of 
960 AU. 




Fig. 4. Representations of the number column density 
of the cloud core derived from simultaneously fitting 
ISOCAM maps at wavelengths of 7 and 15 jum. 



of the x 2 -function. A more detailed error analysis is given 
in Sect. 5. 
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Fig. 5. ISOCAM images of the dense molecular cloud core 
p Oph D at A = 7 /zm (upper left) and A = 15 /jm (lower 
left) and maps resulting from fitting all maps with the 
same number density distribution (right). The white line 
in the upper left panel indicates a spatial length scale of 
about 0.13 arcmin or 0.1 pc (20000 AU). 

4. Determination of constraints on the 3D density 
and temperature distribution 

4.1. The T T -method 

In this section, we want to use a mm map to further 
constrain the density and temperature distribution of the 
dust. However, the dust particles absorbing the radiation 
at 7 and 15 /zm are not obviously the ones dominating 
emission in the mm range. We assume in first approxi- 
mation that the absorption and emission is dominated by 
particles of size a. The absorption in the MIR is described 
by the optical depth along the line of sight 



where B is the Planck function and T is the temperature 
of the dust particles (a detailed description of the mm 
emission is given in Sect. 4.2). The temperature of the 
dust particles can vary with particle size when the dust 
is exposed to strong UV and optical radiation. Pre-stellar 
cores like p Oph D are shielded from the short-wavelength 
part of the spectrum which reduces differences in the dust 
temperature of different dust sizes. For this application, 
we assume a temperature depending on the location but 
not on grain size. Therefore, the emission in the mm is also 
proportional to the optical depth and it is reasonable to 
use the spatial dust particle distribution found by fitting 
maps in the MIR to model the mm emission. 

To analyze the mm map, an important param- 
eter is the temperature of the dust in the core 
with a column density determined from the MIR 
fit. Temperature distributions of embedded and non- 
embedde d spherically symmetric cores ha ve been calcu- 
lated by IStamatellos fc Whitworthl l|2003h ranging from 
6.5 to 18 K. IZucconi et alJ l|200ll) have discussed both ID 
temperature distributions and temperatures in disk-like 
structures. 

The core p Oph D clearly has a more complex struc- 
ture. The temperature of the dust at a certain location x 
within the core is determined by the local power density 
balance of incoming and outgoing radiation 



dX 4tt a 2 Q a6s (A) n(x) B x [T{x)} = 



2tt 



dX it a 2 Qabs(X) n(x) 



1 

4tt 



dcj) I d0sm6 I(\,x,$4>) 



n(x, y, z) dx 



a M c (y,z) = 3 Q 
m d 4p d a 



with the total column dust mass M c , the dust particle 
mass rod, an d the absorption factor Q. For pre-stellar 
cores, the dust sizes are presumably small enough so that 
the Rayleigh limit 2ira <C A is fulfilled for wavelengths 
around 1 mm and commonly assumed chemical composi- 
tion of astrophysical dust. Hence, Q cx a, and the optical 
depth is independent of the dust size as long as the total 
dust mass is not changed. 

The intensity of the mm emission is 

I mm (y,z) oca / n(x,y,z) B\[T(x,y, z)] dx 



with the intensity / and the direction described in spheri- 
cal coordinates 9 and </>. Generally, the radiation intensity 
within the core depends both on the external irradiation 
and the re-emission of the dust particles, coupling / and 
T over all wavelengths, locations, and directions. For a 
complex density structure like the core p Oph D, a 3D 
CRT calculation with self-consistent temperature itera- 
tion is required. However, the numerical effort of these 
codes strongly prohibits any fitting of multi-parameter 
distributions, as the temperature distribution has to be 
Mc(y, z)(Q|ctcrmincd many thousand times in order to find the op- 
timal multi-parameter set in the solution space. 

To estimate the range of the temperatures within 
p Oph D, we calculated the radiation field emerging 
from several density distributions with Gaussians being 
stretched and translated along the line of sight, but with 
a fixed extent and column density in the plane of sky 
to fit the MIR map. The calculations have been per- 
forme d by a 3D CRT code presented inlsteinacker et alJ 
lEioi Csee also ISteinacker et all l2002alb|) . We followed 
IZucconi et all (|200l|) in assuming an embedding envelope 
around the core distribution mimicking the true molecular 
cloud around p Oph D. This outer matter shields the core 
(10) against UV and parts of the optical outer radiation leading 
to one magnitude of optical extinction. As outer radiation 
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field, the mean interstellar radiation field given in Black 
<ll994h. an MIR pow er law distribution as suggested in 
Zuc coni et al . (2001). and the radiation field of the nearby 
B2V star were used IlLiseau et al 1 ll999h . Numerically, the 
direction integral was performed using — 49 equally 
distributed nodes (9j , <pj ) on the u nit sphere along with 
optim ized integration weights ujj l|Steinacker fc Thamml 
Il99fil) . 

Fig. shows results of a first 3D run where we de- 
fined for all Gaussians the HWHM along the x-axis (line 
of sight) to be the mean of the HWHM along the other two 
axis. A priori, there is no reason to believe that the distri- 
bution along x is substantially different from the other two 
axis, so this is an appropriate first guess. The Gaussians 
were not translated along the line of sight in this run 
(ti = 0). The top plot in Fig.Elgfves the density as a func- 
tion of the distance to the southern density peak. As the 
slope of the structure has different gradients in different di- 
rections, the distribution broadens quickly towards larger 
radii, with other maxima at distances of about 20000 AU 
from the southern peak. The optical depths for all grid 
points leading from outside the core to the point x along 
a direction described by the unit vector eg^ and the pa- 
rameterization s is 



70^(0, oo, A) = ct(A) J ds n(x— seg^) . 



(12) 



For A = 7 /im, the resulting dependency on the distance 
of x from the southern density maximum is shown in the 
middle panel of Fig. for all grid cells. Due to the large 
number of cells and directions, we have plotted, for each 
grid cell, only the minimal and maximal optical depth 
when varying the direction. Absorption leads to a max- 
imal optical depth < 2 for photons reaching or crossing 
the densest part of the core, but again the distribution 
of possible optical depth is broad. The bottom panel of 
Fig. El gives the temperature at all grid cells as a func- 
tion of the distance to the density maximum. In the outer 
parts, T reaches some 15 K, while at the maximal den- 
sity, it drops below 9 K. Inbetween, temperatures vary 
strongly depending on the shielding of the different grid 
cells, and no simple approximation of the temperature can 
be derived from this plot. 

In order to perform a fit of the mm emission of the core 
density distribution to the observed mm map, an approxi- 
mative way to find the temperature with moderate numer- 
ical effort is required, but not obvious in view of the large 
spread in the temperature distributions discussed above. 
In the following, we introduce an approximate method to 
determine temperatures for pre-stellar cores. Due to the 
low dust temperatures, the core emits only in the FIR and 
mm. As the core is optically thin to this emission, we ne- 
glect heating of the inner parts by the outer low-density 
parts. The intensity of the illuminating radiation from the 
optical to the FIR, therefore, can be written as 
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Fig. 6. Results from a 3D CRT calculation of the density 
fitting the MIR maps. The third unknown extent of the 
single Gaussian was assumed to be the average of the ex- 
tent in the two other directions. The top panel gives the 
density of the grid cells as a function of the distance from 
the southern density peak. The northern density peaks are 
visible at a distance of about 20000 AU. The middle and 
lower panel show the limits of the optical depth at 7 /im 
and the temperature of the grid cells as function of the 
distance from the southern density peak, respectively. 



with the background intensity Ib g (X). The direction inte- 
gral on the right-hand side of Eq. I|ll|) can be interpreted 
as a weighted averaging over different optical depths, so 
we define a mean optical depth 



N d 



I{X,x, 



I 69 (A)exp [-T0^ (0,oo, A)] 



(13) 



= (A) 



ln { ^2 w i exp [ _Te ^ (°> 00 ' A )] 



(14) 
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For a point in space being surrounded by matter with 
equal column density in all directions, the mean optical 
depth is identical to the optical depth in any direction as 
required. For an embedded point in a core with a gap in 
a certain direction h, we correctly find r(A) ~ T o h .<p h W- 

Fig. gives the relation between the temperature and 
the mean optical depth at 7 /im. Hereinafter, we will call 
this graph the "T— diagram". The relation shows surpris- 
ingly little spread. For low optical depth r < 0.04, the 
dust is exposed to the NIR background radiation field and 
the derived temperature of 15.5 K becomes independent 
of the absorption (dashed line). For larger optical depth 
r > 0.1, the distribution roughly follows a r ~ T 4 -law 
(solid line). This law depends on both the background 
field ib fl (A) and the assumed absorption coefficient k(A). 
In our case, the re-emission of the dust with a tempera- 
ture of 7-15 K will occur in the FIR where the assumed 
absorption law approximately follows k oc A~ 2 . The in- 
tegral on the left-hand side of Eq. I|11J) (representing the 
power density of the radiation emitted by the dust) has 
then an analytical solution oc £(6) T(6) T 6 with the zeta 
function £ and the gamma function T. The power den- 
sity of the incoming radiation on the right-hand side of 
Eq. Ulljl contains the sum of two exponentially-dropping 
contributions from the NIR and FIR peaks of the inter- 
stellar radiation field, and is ~ t 7 ®^ for the background 

field discussed here. This yields the relation T ~ f\'^ m as 
seen in Fig.[7| The T(r)-correlation can be found for other 
wavelengths as well. However, for A < 7 /im, due to the 
high optical depth, the radiation is no longer contributing 
substantially to the heating which broadens the distribu- 
tion. For A > 7 /xm, the relation holds but has more scatter 
as the dust is absorbing less efficiently. 

We will base the following analysis on assuming that 
the relation established in the T— diagram will hold even 
when the detailed 3D structure of the core is changed. The 
key point is that the heating is mainly related to the mean 
optical depth and not so much to the actual spatial dust 
distribution that is leading to this mean optical depth. 
We have tested this assumption by moving and deforming 
the 30 Gaussian 3D distributions representing our core. 
The outcome of the calculations is shown for four typi- 
cal examples in Fig. |SJ To test if the T^relation holds we 
have plotted the HWHM of the distribution of grid cells 
having a temperature T, corresponding to vertical cuts 
through the T^diagram, for a mean optical depth at 7 /jm 
of 0.03, 0.05, 0.07, 0.09, respectively. For convenience, the 
peaks were normalized to 1. The top left panel gives the 
HWHM for the non-translated and non-perturbed den- 
sity distribution as derived from the fit of the MIR map 
and a Gaussian HWHM along the line of sight that is 
the mean of the two other HWHM of the Gaussians. The 
right top panel shows the results for a density distribution 
where we have moved the Gaussians in alternating direc- 
tion along the line of sight. In the lower left panel, the 
Gaussians have been stretched by a factor of 2, and in the 
lower right panel, they were stretched and shifted simul- 
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Fig. 7. Temperature of the grid cells as a function of the 
mean optical depth at 7 ^m. The solid line indicates the 
slope of a power law of the form rT 4 =const. The dashed 
line indicates the maximal temperature of grid cells out- 
side the core but inside the cloud. 

taneously. Evidently, for all four cases, the HWHM of the 
grid cell distribution are around or less than 1 K, with the 
trend that stretching of the Gaussians is decreasing the 
HWHM of the grid cell distribution. The peaks are moved 
by less than 0.6 K compared to the first case. We con- 
clude that the T~ diagram can be used to determine the 
temperature approximately from the mean optical depth 
for all configurations where the Gaussians are shifted and 
stretched within these limiting cases. 

4.2. Fit of the mm map 

To fit the 1.3 mm map of p Oph D, we applied the back- 
ground determination method described already for the 7 
and 15 /zm maps, and obtained a weak background emis- 
sion intensity hg(y, z) in the plane of the sky. The solution 
of the radiative transfer equation including absorption and 
emission along a ray parallel to the x-axis is 

I(+oo,y,z) = 

I bg (-oo,y,z) e -"W"(i/,*) + 

+ 00 

tr(A) J olx' n(x',y,z) B X [T {x 1 , y , z)] e -r,(x',+co, y ,z^ 

— oo 

For wavelengths in the mm range, the optical depth will be 
small (t x < 10 ) and the exponentials are approximately 
1. The solution Eq. I|15[) contains the unknown number 
density along the ray n as well as the temperature distri- 
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Fig. 9. Original 1.3 mm map (left), theoretical map for a column density obtained from the MIR-map fit and a 
constant temperature of T=10 K (middle), and the fitted map using the T— relation (right). The grey circle in the left 
image represents the beam size. 
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Fig. 8. Normalized grid cell number at 7 /Ltm as a func- 
tion of the temperature for an optical depth of 0.03, 0.05, 
0.07, and 0.09, respectively, for 4 different configurations 
sketched in the inset (30 unperturbed Gaussians, shifted, 
stretched, and shifted/stretched). The HWHM can be 
read from the dashed line. 



bution T. Assuming 3D Gaussians with 30 main axes and 
30 relative translations along the line of sight, we know 
the full 3D density n. The mean optical depth T7 Mm can 



be found and translated into a temperature using the T— 
diagram by ray-tracing from different directions. Finally, 
we integrate along the line of sight according to Eq. i|14[l to 
find the intensity. By varying the translation and HWHM 
of the Gaussians along the line of sight with the simulated 
annealing technique, we perform a fit of the mm map. 

The resulting fit to the mm map is shown in Fig. |5J 
The left image shows the observed map. Assuming a con- 
stant dust temperature of 10 K, the column density map 
derived from the MIR image can be used to calculate a 
corresponding mm map (middle image). Aside from the 
noise of the map, comparison of the map reveals several 
differences. While the observed map has a more bar-like 
structure with an approximately constant flux level along 
the middle of the bar, the mm map for constant tem- 
perature has a clear gap between the two main emission 
maxima seen in MIR. Moreover, the peak flux value of 
the northern maximum is lower compared to the maxi- 
mum flux in the lower part of the image. Allowing the 
program to vary the extent of the Gaussians along the 
line of sight, the mm emission can be increased by stretch- 
ing individual clumps. This way, the optical depth to the 
considered point in space is lowered allowing a more ef- 
ficient heating and a rise in temperature. While the den- 
sity is slightly lowered by the stretching, the overall effect 
is to increase the emissivity. The optimization program 
therefore stretched the Gaussians in regions where the ob- 
served image was brighter, especially in-between the two 
MIR-emission maxima. Also, the Gaussians forming the 
maximum in the northern part of the image have been 
stretched and moved away from each other to decrease 
the mean optical depth and to increase the emissivity. 
The right image in Fig. El shows the results of shifting 
and stretching/compressing the Gaussians along the line 
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of sight to fit the mm map. The emission pattern now is 
more bar-like, but still a gap is visible that cannot be seen 
in the observed picture. The brightness of the two maxima 
is equalized, and flux is moved from the southern emission 
maximum towards the middle structure. 

Discussing the deviations, one important point has to 
be raised in the beginning. A much better fit could be 
achieved when the Gaussians would be allowed to move in 
the plane of sky. But as the method is designed to fit both 
the MIR and mm maps, the Gaussians are fixed in the 
plane of the sky to enforce a simultaneous fit of the MIR 
map. This includes that the total mass of each Gaussian 
is fixed. 

The freedom in the fit is reduced due to the limited 
number of Gaussians, and a better fitting especially in- 
between the two MIR-emission maxima could be expected 
from a run with more Gaussians. 3D radiative transfer, 
even when using the T^method and combined use of sim- 
ulated annealing is at the limit of the currently avail- 
able computer powers so that a substantial increase of 
the number of Gaussians is not feasible. In addition, the 
low-resolution mm map used for the fit has a sub-structure 
that prohibits unambiguous identification of fine-structure 
in the image. We hope to improve this with upcoming 
new mm maps. A deeper discussion of the various uncer- 
tainties, errors, and approximations is given in the next 
section. 

In Fig. we have visualized the 3D dust density 
structure of the core by showing two iso-density surfaces 
at 5% (solid surface) and 1 % (semi-transparent surface) 
of the maximal density and rotating it against an artifical 
background. The lower density maximum is compact with 
an extension to model the MIR emission seen in the left 
lower part of the MIR image. The northern structure is 
more complex and extended to model both the broader 
structure seen in the MIR and the high emissivity seen in 
the mm. In Tab^ we summarize the overall properties of 
the cloud. 

The 3D dust temperature distribution is visualized in 
Fig.^Jby showing iso-temperature surfaces for increasing 
temperature. The left panel gives the view along the line 
of sight (corresponding to the top left image in Fig. I1U|) , 

Table 1. Derived overall properties of the core p Oph D 



Spatial extent 0.22 pc 

Total gas mass 2.3 Mq 

Total gas mass of northern 0.8 Mq 
peak 

Total gas mass of southern 1.5 Mq 
peak 

Maximal density at southern 30 m -3 
peak 

Main axes aspect ratio 2:2:1 

Mean gas density 8 x 10 -16 kg m 

Temperature range 10 — 15 K 




Fig. 10. Iso-density surfaces of the 3D dust density struc- 
ture of the cloud core p Oph D at 5% (solid) and 1 % (semi- 
transparent) of the maximal density for different viewing 
angles ranging from to 320 degrees with respect to the 
line of sight seen against an artificial background. 

while in the right panel, the configuration is rotated by 
90° in a plane perpendicular to the plane of the sky 
(corresponding to the image with the label "135" in 
Fig. 110(1 . The iso-temperatures are 12 K, 12.8 K, and 
13.8 K, respectively (top to bottom). In general, the dis- 
tribution follows the density distribution, as the densest 
parts are well-shielded and have low temperatures. In 
the northern part, however, there is a low-temperature 
region that is not only caused by shielding due to 
matter around the local density maximum, but also by 
other nearby local density maxima. This effect becomes 
important when the density distribution consists of many 
clumps with comparable maximal density values as in the 
northern clump. We note that such a realistic tempera- 
ture distribution including shielding effects can only be 



12 



Steinacker et al: 3D Structure of p Oph D 











Fig. 11. Iso-tcmpcraturc surfaces of the 3D dust temper- 
ature structure of the cloud core p Oph D. The left panel 
shows the distribution seen along the line of sight, while, 
for the right panel, the view has been rotated by 90°. The 
iso-temperatures are 12 K, 12.8 K, and 13.8 K, respec- 
tively (top to bottom). 



obtained doing RT calculations, but that the proposed 
fast T~ method includes this effect. Animations show- 
ing the density and temperature structure are displayed at 
http: / / www.mpia-hd.mpg.de /homes/stein/ Ani/movief. htm 



5. Ambiguity and error discussion for the 3D 
modeling 

Deriving the complex 3D structure of pre-stellar cores 
from projected images is a difficult and challenging task. 
The method and results presented in this paper are a 
first step to enter this rich 3D world. While the advan- 
tage over a ID or 2D approach is obvious in view of the 
complex structures of observed cloud cores, the various 
errors arising from approximations and limited knowledge 
of the physical conditions make a careful error analysis in- 
evitable. Lacking any analytical solution, however, a quan- 
titative error estimate is often not possible. In turn, many 
of these errors or ambiguities might be reduced when maps 
of higher resolution, sensitivity, and dynamical range at 
the considered wavelengths or images at other wavelengths 
become available (e.g. by SPITZER, HERSCHEL, JWST, 
ALMA), or when independent information, e.g., about the 
background radiation field can be added to the currently 
available information. 



In the following, we will discuss the ambiguities and er- 
rors arising from usage of the numerical algorithms, from 
the approximations and assumptions, and from the obser- 
vational errors. 



5.1. Numerical errors 

The column density and emission integrations for given 
density and temperature distributions are performed with 
a 5th-order Runge-Kutta scheme incorporating adaptive 
step size control. The error of this ray-tracing can be ne- 
glected compared to the other sources of error. 

The density distribution contains 30 three-dimensional 
Gaussians. In Sect. 6, we already concluded that the num- 
ber of Gaussians might be too low to resolve the material 
between the two MIR peaks sufficiently, causing deviations 
in the mm emission from this region of the core. We did 
extensive tests including Gaussians with arbitrary orien- 
tation in space. Indeed, a smaller number of Gaussians of 
this type is required to fit a density distribution than for 
Gaussians with axes along the Cartesian axes. However, 
the formula for Gaussians with arbitrary orientation in 
space is substantially more complex increasing the run 
time of the fitting routine beyond current computer ca- 
pabilities. We are currently elaborating, if the code can 
be accelerated by other means to be able to increase the 
number of Gaussians. 

The 3D CRT code has been tested in the framework 
of the 2D spectral benchmark published in lPascucci et alJ 
(2004|), the first benchmark to test five Monte-Carlo and 
grid/ray-tracing CRT codes by calculating the spectral en- 
ergy distribution of a star surrounded by a circumstellar 
dust disk. The deviations of the results for optical depths 
considered here were below 10%, and most of the devia- 
tions were caused by the scattering term, whereas scatter- 
ing plays no important role for the heating of pre-stellar 
cores. Nevertheless, we point out that there is no 2D CRT 
benchmark comparing images nor a 3D CRT benchmark 
of any kind. 

The limitations and possible errors of the proposed 
2V method have been discussed already in Sect. 4 for the 
application to the cloud core p Oph D with an approx- 
imate total error of about 1 K or less. In general, the 
method requires to rerun the full 3D CRT code for any 
newly considered background radiation field in order to 
obtain the T— relation. In cases of strong sub-structure 
of the core material, the number of considered directions 
for the averaged optical depth might have to be increased. 
The possible errors arising from additional, but unresolved 
shielding can be neglected in most cases, but for structures 
with holes, the method might underestimate the tempera- 
ture when missing to resolve the illumination through the 
hole. Comparing the scale of the sub-clumps with the spa- 
tial resolution of 1000 AU/pixel for cores in p Oph with 
MIR-instruments like ISOCAM, it is evident, though, that 
such a strong sub-clumping is currently not observed. 
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Operating in a complex 210-parameter space, the main 
question is the ambiguity of the derived optimized pa- 
rameter set. Optimization algorithms like the gradient 
search usually succeed in finding a local minimum, but 
often fail in localizing the deepest minimum of the y 2 - 
function of the fit. For this purpose, more sophisticated 
methods like simulated annealing or genetic algorithms 
are applied which have a higher probability in finding the 
deepest minimum. A prove for completeness is impossible 
in our ca se, but from our former experience in applying the 
method llThamm et al1ll994t ISteinacker fe Thamml ll996: 
ISteinackeTi§r^eTmhiell2003j) as well as from exploring dif- 
ferent initial parameter sets, we are confident that the so- 
lutions are close to the density structure presented here. 
Nevertheless, it is correct to state that we can only show 
that the derived density structure simultaneously fits the 
three images discussed here, but that we can not exclude 
that there is another distribution with a comparable \ 2 - 
value. We note that the presented density distribution is 
a valuable working hypothesis for further studies and that 
we intend to remove some of the ambiguity when more 
data become available. 

The proposed method has another ambiguity that is 
common in interpreting projected images of 3D configu- 
rations. As soon as Gaussians are separated along the line 
of sight (e.g. when the distance of the maxima is much 
larger than the HWHM) and not indirectly connected by 
other Gaussian, there is no way to determine the distance 
of the clumps. This is because the mm fit uses the temper- 
ature difference arising from overlapping Gaussians, and 
therefore, the temperature is constant for all clumps of 
large separation. Moreover, separated clumps can be ex- 
changed in their position along the line of sight without 
changing the observed MIR and mm flux. 

For the density structure derived here, this is the case 
only for a few low-density Gaussians. Hence, we conclude 
that this ambiguity does not affect the derived distribution 
of p Oph D. 

5.2. Approximation and assumption errors 

The assumptions concerning dust properties have already 
been discussed in Sect. 3 and the beginning of Sect. 4. 
As the total dust mass derived from the MIR maps does 
not strongly depend on the particle size, we expect the 
uncertainties in the optical properties of the dust and the 
gas/dust ratio to dominate the error in the mass. It should 
be noted that another small error has been introduced 
by comparing a flux that has been observed through the 
wavelength filters of ISOCAM to images calculated for 
a single wavelength. An advanced modeling of the dust 
properties will be subject to a forthcoming project. 

It might be objected that the expansion of the den- 
sity distribution may cause substantial errors when fit- 
ting filamentary structures, especially when inclined to the 
Cartesian axes. Firstly, the images with the 1000 AU/pixel 
resolution do not resolve such strongly elongated struc- 



tures. Secondly, in order to avoid un-physically flat struc- 
tures, the fitting algorithm allows main axes ratios of up 
to 10. The maximal ratio obtained in the fits is about 5. 
Moreover, elongated structures inclined to the Cartesian 
axes are fitted by several Gaussians with moderate main 
axes ratios. We therefore argue that the expansion into 
3D Gaussians does not produce large fitting errors. 

For the determination of the column density from the 
MIR images, we assume that the invisible background be- 
hind the core can be extrapolated from the visible back- 
ground in the vicinity of the core (the details of the used 
2 step-interpolation are given in Sect. 3). The stability 
of the method was tested by varying interpolation radius 
and core area. The local variations within the interpo- 
lation radius are of the order of 10-15% arising either 
from variations in the illuminating MIR flux of the photo- 
dissociation region or from low-density molecular cloud 
structures along the line of sight. 

The exact temperature distribution depends on the as- 
sumed properties of the external illuminating radiation 
field. The 3D radiation field for a core within p Oph is not 
known precisely, as the position within the cloud along the 
line of sight is unknown, and the MIR flux of the photo- 
dissociation region and nearby other cores has not been 
modeled in detail. We note that our approach to have a 
composition of interstellar radiation field, MIR power law 
component and nearby B2V star is more realistic than a 
ID background field as it is assumed in ID fits. 

5.3. Observational errors 

The noise of the MIR images is about 0.3 MJy/sr for the 
LW2 ISOCAM filter, so that the MIR flux is uncertain 
to a relative error of less than 10 %. For the mm image, 
the noise is 6 mJy/13" beam for the core center and the 
relative flux error is of the order of 15 %. 

The simultaneous fit of different images requires an ac- 
curate pointing. This is especially valid for the comparison 
of the MIR and mm image, as the Gaussians are not free 
to move in the plane of sky when performing the fit to 
the mm map. An offset would cause the fitting procedure 
to compensate deviations in position by varying the tem- 
perature un-physically. The absolute pointing error in the 
images was less than 1.4" (2a) for the ISO telescope and 
3" (lcr) for the IRAM 30m telescope, which is less than 
the typically obtained HWHM of the Gaussians. 

The technique used for the IRAM 30m 1.3 mm bolome- 
ter observations ("dual beam") is known to filter ou t 
some of the ext e nded emissions (e.s. lAndre et al.l lT996"). 
iMotte fc Andrei l|200lh have simulated the effects of the 
dual-beam technique on model objects with intensity pro- 
files following simple power laws and calculated the loss of 
flux arising in these observations. The flux loss is greater 
for low power-law indices and increases with increasing 
radius. We used their modeling to evaluate the flux loss 
in p Oph D and found it to be negligible. 
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In Sect. 3, we have used a zodiacal light component 
being constant over the observed area which has an un- 
certainty of about 20%. Due to the large contribution of 
the zodiacal light to the total flux and its large error, 
this foreground flux is dominating the uncertainties in the 
modeling that arise from observational errors. 

6. Comparison with 3D density distributions from 
gravo-turbulent star-formation models 

In the picture of gravo-turbulent star formation, molec- 
ular cloud cores are not quasi-static equilibrium objects. 
Rather they are dynamically evolving as part of the over- 
all turbulent cascade. The same complex and stochastic 
flow that forms a core at first place continuously reshapes 
its structure and even may disperse it again. Those cores 
with an excess of gravitational energy collapse rapidly to 
form stars, while the others with sufficiently large internal 
or kinetic ene rgies re-expand once the turbulent compres- 
sion s ubsides l|Tavlor et, alll99fil:IVazauez-Semadeni et all 
120041) . Consequently, molecular cloud cores are observed 
having a variety of different shapes, ranging from highly 
elongated filaments via cometary-shaped struct ures to 
very roundish spherically symmetric cores (see e.g. iMversI 
1999). This sequence roughly reflects the increasing im- 
portance of self-gravity compared to turbulent kinetic en- 
ergy. Objects that are dominated by turbulent ram pres- 
sure tend to have on average more complex morphological 
appearance than more quiescent cores on the verge of grav- 
itational collapse l|Klessen et alJliooih . The latter usually 
have a rather generic density structure with flat inner pro- 
file, followed by an approximate power-law fall-off further 
out, and apparent truncation at some maximum radius. 

The double-peaked and elongated nature of the core 
p Oph D places it somewhere in between these two ex- 
tremes. Its 3D density structure is typical for pre-stellar 
cores in their early phases of evolution when the dynamics 
is still domi nated by external co mpression rather than by 
self-gravity l|Klessen et al.l2004l) . The relatively small den- 
sity contrast of about 30 is additional evidence for that. 
To illustrate this point, Fig. E| shows the column den- 
sity map of a typical pre-stellar core in a smooth parti- 
cle hydrodynamics (SPH) simulation of a self-gravitating 
supersonically turbulent molecular cloud for comparison 
(integration along the Cartesian axes). 

The computation focuses on a cubic volume within 
the cloud and follows the evolu tion using smooth par- 
ticle hydrodynamics (SPH; see iBend 1199ft iMonaghanl 
119921) . This is a particle-based scheme to solve the equa- 
tions of hydrodynamics well-suited for complex flows 
with large density contrasts. Turbulence is continuously 
driven on large scales using a non-local Gaussian field 
such t hat the root m ean square Mach number is con- 
stant llMac Low! [1999) with value M. sa 6. Details of 
the numerical method and associated reso l ution issues 
are discussed , e.g. . by iKlessen k, BurkerU l)2000j) and 
iKlessen et alJ l|2000l) . The gas is taken to be isothermal 
with temperature T = 11.3 K, corresponding to a sound 



speed c s — 0.2 km s -1 . The computational domain corre- 
sponds to a cube of 2.5 pc length, containing 1140 M® of 
molecular cloud gas. Note, Figure 12 shows only a fraction 
of 1/103 of the total volume centered on the considered 
protostellar core. Compl ementary aspects of the ad opted 
model are discussed b v iBalle steros-Paredes e t alJ (their 
model LSD at time U .120031) .ISchmeia fc Klessenl (m odel 
M6k2.l2004l.lJaPDsen fc Klessenl (model M6k2. l200^ . and 
IKlessen et al I Cmodel LSD. I2004I) . 

The model core in Fig. is selected because its mor- 
phological appearance closely resembles that of p Oph D. 
At the depicted time it is gravitationally unstable and 
will collapse to build up a protostar during its subsequent 
evolution. However, due to the highly stochastic nature 
of compressible turbulent flows the fate of a molecular 
cloud core is predictable only in a statistical sense. The 
selected model core, therefore, can at best be indicative 
of what may happen to p Oph D in the next 10 5 years. 
Nevertheless, it is fair to speculate that gas around the 
southern density maximum will go into collapse. The fate 
of the more extended and less dense northern core is hard 
to predict without further study of possible infall motion. 
When collapsing into two fragments, the core p Oph D 
would give rise to a wide binary system. A detailed inves- 
tigation of the dynamical evolution of the core based on 
the observational data will be the subject of future studies. 
This will require a combined modeling of the continuum 
and line data to incorporate velocity informations. 

7. Conclusion and Outlook 

In this paper, we propose a new method to model the 
three-dimensional dust density and temperature structure 
of dense molecular cloud cores. The method is based on 
the fits of multiple continuum images in which the core 
is seen in absorption and emission and requires only a 
few computations with a 3D CRT code. The large pa- 
rameter space of a three-dimensional structure is covered 
using simulated annealing as optimization algorithm. One 
of the key points of the proposed method is the use of a 
T— relation that has been derived from the 3D CRT runs. 
It links the mean optical depth at a given wavelength from 
outside to a spatial point within the core to the dust tem- 
perature at this point, making a fast evaluation of the mm 
emission along the line of sight possible. 

We have applied the method to model the dense molec- 
ular cloud core p Oph D. In the MIR, the core is seen in 
absorption against a bright background from the photo- 
dissociation region of the nearby B2V star, illuminating 
the cloud from behind. Two ISOCAM images at 7 and 15 
jLtm have been fitted simultaneously by representing the 
dust distribution in the core with a series of 3D Gaussian 
density profiles. The background emission behind the core 
was interpolated from nearby regions of low extinction. 
Using simulated annealing, we have obtained a 2D column 
density map of the core. The column density of the core 
has a complex elongated pattern with two peaks, with the 
southern peak being more compact. To retrieve the full 
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Fig. 12. Column density map of a typical pre-stellar core based on a SPH simulation of a self-gravitating supersonically 
turbulent molecular cloud. The density is integrated along the three Cartesian axes. 



three-dimensional structural information, we have calcu- 
lated the temperature structure of the core with a 3D CRT 
code, assuming several limiting cases of core extent along 
the line of sight. It was shown that, for a given position in 
the core, the relation between the dust temperature and 
the mean optical depth from outside to this position varies 
within less than 1 K, when changing the shape of the core 
along the line of sight. Therefore, the T^relation was used 
to estimate quickly the temperature at a given point by a 
fast direction averaging over the optical depth at just one 
wavelength. Using this mean, we have varied extent and 
position of the Gaussian components of the density along 
the line of sight until a fit of both the MIR and a 1.3mm 
map, obtained with the IRAM 30m telescope, could be 
achieved. We have presented the three-dimensional dust 
density and temperature structure, revealing a condensed 
southern part and a more extended and complex northern 
part. We have addressed several sources of errors, namely 
the background determination method, the assumed dust 
particle properties, errors in the maps, the ambiguity in 
the derived distribution, and the approximation error of 
the Tf method. As we did not reach a perfect fit of the 
mm image between the two maxima due to the low num- 
ber of Gaussians, we expect that the detailed 3D structure 
in-between the density maxima might change when using 
a higher number of Gaussians. The general structure of 
a condensed southern maximum and a complex northern 
multi-clump region, however, is clearly resolved within the 
limits of the presented analysis. This structure is in gen- 
eral agreement with recent gravo-turbulent collapse calcu- 
lations for molecular clouds. We speculate that the south- 
ern condensation belongs to the category of cloud clumps 
which are dominated by self-gravity. It may be collapsing 
to form a star in the near future. The northern clumps, 
however, may still have sufficiently large internal or kinetic 
energy to re-expand and merge into the ambient medium 
once the turbulent compression subsides. 

In this application, we determine about 210 free pa- 
rameter to fit about 2000 pixels simultaneously, but hid- 



den in the assumptions are many more free parameter 
like the dust properties or details of the illuminating back- 
ground. It has to be kept in mind, though, that the hidden 
parameters are present in any other, more simple modeling 
as well, and often with more unrealistic approximations. 
The p Oph cloud is illuminated from behind in contradic- 
tion to a ID boundary condition for the incoming radia- 
tion. 

The perspective of applying the method is promising. 
There are a number of well-observed cores where images 
at more than two wavelengths are available. With every 
additional image, the ratio of new unknown background 
parameter versus constrained parameter decreases, as the 
density parameters are unchanged. This will increase the 
accuracy of the determined density structure. Moreover, 
the derived TV-relations for p Oph D can be used in other 
applications incorporating the temperature. 

The ultimate goal of applying the method to well- 
observed cores, however, will be to address the key ques- 
tion of early star formation, namely if the considered cores 
have in-falling material. The current line observations pro- 
vide the molecular line emission flux integrated over all 
moving gas cells along the line of sight. In the general case, 
gas motion and emissivity of the cells can not be disen- 
tangled, and the ID approximation or shearing layers are 
assumed to unfold it. Without unfolding it, infall motion 
can be mixed up with rotational motion leaving it unde- 
cided if the core shows any sign for the onset of star for- 
mation. This is changed if the core has been investigated 
with the TV-method. Knowing the full 3D structure in dust 
density and temperature, the line of sight-integral can be 
inverted providing the complete kincmatical information 
if the considered line is optically thin and a model for the 
depletion of the considered molecules is used. This direct 
verification of infall motion will be subject to a forthcom- 
ing publication. 
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